### File    : replication_empirics.R
##  Author  : Erin York
##  Notes   : Replication materials for Slough, York, and Ting -
#             A Dynamic Model of Primaries - empirics
#             R version 3.5.2; Loaded package versions: gridExtra_2.3, bindrcpp_0.2.2,
#             RColorBrewer_1.1-2, ggplot2_3.1.0, xtable_1.8-3, stargazer_5.2.2, dplyr_0.7.8 
#             plyr_1.8.4, beepr_1.3     

rm(list = ls())

library(plyr)
library(dplyr)
library(beepr)
library(stargazer)
library(xtable)
library(ggplot2)
library(RColorBrewer)
library(gridExtra)

load("05_replication_data.Rdata")



# Merge datasets ----------------------------------------------------------

# trim DALP variables
parties<- dalp %>% select(1:6, d1, d2, d3, d4, d5, dw, starts_with("cosalpo"), partysize)

# filter latam data to relevant years & add DALP measures
polar<- latam %>%
  filter(year >= 2006, year <= 2010) %>%
  left_join(parties, by = c("country", "partyDALP" = "party"))

# n parties per election
polar<- polar %>%
  group_by(country, year) %>%
  mutate(party_dat = sum(!is.na(d1))) %>%
  ungroup()

# create polarization metrics
# polar 1 = ideological variance
# polar 2 = ideological range
polar<- filter(polar, party_dat > 1) %>%
  group_by(country, year) %>%
  mutate(polar1 = sd(dw, na.rm = T),
         polar2 = max(dw, na.rm = T) - min(dw, na.rm = T)) %>%
  ungroup() %>%
  mutate(primary = as.numeric(primary) - 1) %>%
  filter(!is.na(dw))

polar<- arrange(polar, polar1) %>%
  mutate(Primary = mapvalues(as.factor(primary), from = c(0,1), to = c("No", "Yes")))

polar2<- polar %>%
  dplyr::group_by(Primary) %>%
  summarise(polarization = mean(polar1),
            marginofwin = mean(margin))

# average polarization and vote margine w/ and w/o primaries
polarprim<- mean(polar$polar2[polar$primary == 1])
polarnoprim<- mean(polar$polar2[polar$primary == 0])
marginprim<- mean(polar$margin[polar$primary == 1])
marginnoprim<- mean(polar$margin[polar$primary == 0])



# Figure 1 ----------------------------------------------------------------
### Primary election laws 

laws$country <- with(laws, factor(country, levels = rev(levels(country))))
laws$Primary_Laws <- with(laws, factor(status, 
                                       levels = c("Unspecified", "Voluntary", 
                                                  "Mandatory, Not Enforced",
                                                  "Mandatory, Enforced")))
laws$Primary_Held <- ifelse(is.na(laws$election.year), NA, as.character(laws$election.year))
k <- which(laws$Primary_Held %in% c("One Party", "Some Major Parties", "Some Parties"))
laws$Primary_Held[k] <- "Some Major Parties"
laws$Primary_Held <- factor(laws$Primary_Held, levels = rev(c("None", "Some Major Parties", "All Major Parties")))

laws_ss <- filter(laws, !is.na(Primary_Held))

a <- ggplot(NULL, aes(x = year, y = country), size = 1) + 
  geom_tile(data = filter(laws, !is.na(Primary_Laws)), aes(fill = Primary_Laws), size = 0) + 
  geom_point(data = laws_ss, aes(shape = Primary_Held), size = 4, bg= "gray85") + theme_bw() + 
  xlab("Year") + ylab("") + ggtitle("Primary Laws in Latin America by Country, Year") + 
  scale_shape_manual("Primaries Held in \nPresidential Election", labels = c("All Major Parties", "Some Major Parties", "None"), values = c(24, 22, 21)) +
  # scale_colour_manual("Primaries Held in \nPresidential Election", labels = c("All Major Parties", "Some Major Parties", "None"), values = c("gray75", "gray50", "gray25")) +
  scale_fill_manual("Primary Laws", values = c("#a1dab4", "#41b6c4", "#2c7fb8", "#225ea8"), na.value = "transparent") +
  theme(plot.background = element_rect(fill = "transparent",colour = NA))
ggsave(a, filename = "Output/Figure_1_primary_laws.pdf", width = 10, height = 5.25)



# Figure 2 ----------------------------------------------------------------
### Primary adoption in presidential elections as a function of polarization


b <- polar %>%
  mutate(prim2 = mapvalues(primary, from = c(0, 1), to = c("No Primary", "Primary"))) %>%
  group_by(prim2) %>%
  mutate(mean = mean(polar2),
         med = median(polar2)) %>%
  ungroup() %>%
  mutate(lab = ifelse(prim2 == "No Primary", "No Primary (n = 30)", "Primary (n = 13)")) %>%
  ggplot(aes(x = polar2)) + geom_histogram(aes(y=..count../tapply(..count..,..PANEL..,sum)[..PANEL..]), bins = 30, col = "white") + 
  facet_grid(lab~.) +
  theme_minimal() +
  geom_vline(aes(xintercept = mean), lwd = 1.2, col = "red", lty = 2) + 
  xlab("Polarization") +
  ylab('Proportion (by panel)') 
b

ggsave(b, filename = "Output/Figure_2_polarization.pdf", width = 6, height = 4)



# Figure 3 ----------------------------------------------------------------
### Party-level patterns of primary adoption in four Latin American countries

primaries <- latam %>% 
  mutate(Party = partyname, 
         Country = country,
         Primary_Held = primary,
         Vote_Share = voteshare)

# Colombia data frame
col <- filter(primaries, Country == 'COLOMBIA') %>% 
  dplyr::group_by(Party) %>% mutate(y1 = -min(year))
col$party1 <- factor(col$partyname) 
levels(col$party1) <- c(levels(col$party1)[1:4], "PC/P de la U", "PC/P de la U", "PV")
col$Party <- factor(col$party1, levels = levels(col$party1)[c(1, 3, 6, 5, 4, 2)])
col$Party <- reorder(col$Party, col$y1)

seg <- group_by(col, year) %>%
  arrange(year, desc(voteshare)) %>% filter(Party != "CD")

seg2 <- aggregate(seg, list(seg$year), FUN = head, 1) %>% 
  select(year, Country, Party, Vote_Share, Primary_Held, y1) %>%
  mutate(xstart = as.numeric(as.character(year)), xend = xstart + 4)
seg2$Party <- factor(seg2$Party, levels = levels(col$Party))

# Brazil data frame
bra <- filter(primaries, Country == 'BRAZIL') %>%
  dplyr::group_by(Party) %>% mutate(y1 = -min(year))
bra$Party <- factor(bra$Party, levels = unique(bra$Party))
bra$Party <- reorder(bra$Party, bra$y1)

segb <- group_by(bra, year) %>%
  arrange(year, desc(voteshare)) 

segb2 <- aggregate(segb, list(segb$year), FUN = head, 1) %>% select(year, Country, Party, Vote_Share, Primary_Held, y1) %>%
  mutate(xstart = as.numeric(as.character(year)), xend = xstart + 4)
segb2[1, 8] <- 1994
segb2[7, 8] <- 2016
segb2$Party <- factor(segb2$Party, levels = levels(bra$Party))

# Honduras data frame
hon <- filter(primaries, Country == 'HONDURAS')%>%
  dplyr::group_by(Party) %>% mutate(y1 = -min(year))
hon$Party <- factor(hon$Party, levels = unique(hon$Party))
hon$Party <- reorder(hon$Party, hon$y1)

segh <- group_by(hon, year) %>%
  arrange(year, desc(voteshare)) 
segh2 <- aggregate(segh, list(segh$year), FUN = head, 1) %>% select(year, Country, Party, Vote_Share, Primary_Held, y1) %>%
  mutate(xstart = as.numeric(as.character(year)), xend = xstart + 4)
segh2b <- data.frame(year = c(2005, 2009, 2013, 2017), Country = "HONDURAS", Party = "PN", Vote_Share=100, Primary_Held = NA, xstart = 2005, xend = 2005, y1 = -1981)
segh2 <- rbind(segh2, segh2b)
segh2$Party <- factor(segh2$Party, levels = levels(hon$Party))

# Costa Rica data frame
cr <- filter(primaries, Country == 'COSTA RICA')%>%
  dplyr::group_by(Party) %>% mutate(y1 = -min(year))
cr$Party <- factor(cr$Party, levels = unique(cr$Party)[c(5, 4, 3, 2, 1)])
cr$Party <- reorder(cr$Party, cr$y1)

segcr <- group_by(cr, year) %>%
  arrange(year, desc(voteshare)) 
segcr2 <- aggregate(segcr, list(segcr$year), FUN = head, 1) %>% select(year, Country, Party, Vote_Share, Primary_Held, y1) %>%
  mutate(xstart = as.numeric(as.character(year)), xend = xstart + 4)
segcr2$Party <- factor(segcr2$Party, levels = levels(cr$Party))

col <- select(col, -party1)

d1 <- rbind(col, bra, hon, cr)
d1$Party <- factor(d1$Party, levels = c(levels(col$Party), levels(bra$Party), levels(cr$Party), levels(hon$Party)[2]))
d1$Primary_Held <- ifelse(d1$Primary_Held == "yes", "Yes", "No")

d2 <- rbind(seg2, segb2, segh2, segcr2)
d2$Party <- factor(d2$Party, levels = levels(d1$Party))


d3 <- filter(d1, Vote_Share > 20)
d3$Primary_Held <- ifelse(d3$Primary_Held == "Yes", "Yes", "No")

p <- ggplot(d3, aes(x = Party, y = year)) + ylim(c(1970, 2018)) + 
  geom_hline(data = d1, aes(yintercept = year), col = "grey", lty = 3) + coord_flip()

q <- geom_segment(data = d2, aes(y = xstart, yend = xend, x = Party, xend = Party), lwd = 2, 
                  col = "grey60", alpha = 0.5) 

c <- p + geom_point(aes(shape = Primary_Held, col = Primary_Held), size = 4) + theme_bw() + 
  scale_color_brewer("Primary Held", palette = "Set2") + scale_shape_discrete("Primary Held") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + facet_grid(Country ~., scales = "free_y") + xlab("") + 
  q + coord_flip() + ggtitle("Primary Adoption and Incumbency in Four Countries, 1974-2014, \nAmong Parties Receiving a Vote Share > 20%")


ggsave(c, filename = "Output/Figure_3_party_adoption.pdf", width = 7, height = 8)


# APPENDIX TABLE 1 --------------------------------------------------------
### Polarization and use of primaries

# models 1 and 2
a2 <- lm(primary ~ polar1 , data = polar)
a21 <- lm(primary ~ polar1 + year + partysize + winvote, data = polar)

# models 3 and 4
b2 <- lm(primary ~ polar2 , data = polar)
b21 <- lm(primary ~ polar2 + year + partysize + winvote, data = polar)


# function to construct bootstrapped standard error
block_bs_polar <- function(){
  nclust <- length(unique(polar$country))
  levels <- unique(polar$country)
  samp <- as.vector(sample(levels, size = nclust, replace = T))
  orig <- polar$country
  dat <- ldply(lapply(samp, FUN = function(x){polar[orig == x,]}), data.frame)
  a2 <- lm(primary ~ polar1 , data = dat)
  a21 <- lm(primary ~ polar1 + year + partysize + winvote, data = dat)
  b2 <- lm(primary ~ polar2 , data = dat)
  b21 <- lm(primary ~ polar2 + year + partysize + winvote , data = dat)
  return(c(coef(a2)[1:2], coef(a21)[1:5], coef(b2)[1:2], coef(b21)[1:5]))
}

set.seed(1234)
bs3 <- replicate(n = 2000, expr = block_bs_polar())

# return bootstrapped standard errors
se <- apply(bs3, MARGIN = 1, FUN = sd)

sea2 <- se[1:2]
sea21 <- se[3:7]
seb2 <- se[8:9]
seb21 <- se[10:14]

mean(polar$polar1)
sd(polar$polar1)
mean(polar$polar2)
sd(polar$polar2)

# print table 2 in latex format
sink("Output/Table_A1_polarization.tex")
stargazer(a2, a21, b2, b21,
          se = list( sea2, sea21, seb2, seb21),
          order = c("polar1", "polar2", "year", "partysize", "winvote", "Constant"),
          dep.var.labels = c( "Probability of Holding a Primary"),
          covariate.labels = c("Ideological Variance", "Ideological Range",
                               "Year", "Party Size", "Winning Voteshare"),
          
          omit = c("country+", "Constant"),
          add.lines = list(c("Prop. of Obs. with Primaries ",rep(".30", 4)),
                           c("Ideological Variance Mean (SD)", rep("2.2 (1.0)", 2), rep("", 2)),
                           c("Ideological Range Mean (SD)", rep("", 2), rep("4.2 (1.8)", 2))),
          omit.stat = c("f", "ser"),
          #  title = "Polarization - Party-Year",
          float = F)
sink()



# APPENDIX TABLE 2 --------------------------------------------------------
### Primaries and Electoral Prospects


latam2<- latam %>% group_by(country, year) %>%
  mutate(nparties = length(unique(position))) %>%
  arrange(country, year, desc(voteshare)) %>%
  mutate(vs1 = voteshare[1],
         vs2 = voteshare[2],
         marg = ifelse(voteshare == vs1, voteshare - vs2, voteshare - vs1)) 

latam2<- mutate(latam2, election = paste0(country, "_", year),
                primary_numeric = as.numeric(as.character(mapvalues(primary, from = c("no", "yes"),
                                                                    to = c(0, 1)))))

vote1<- (lm(voteshare ~ primary, latam2))
vote2<- (lm(voteshare ~ primary + nparties + Decade + country, latam2))


block_bs_voteshare <- function(){
  nclust <- length(unique(latam2$election))
  levels <- unique(latam2$election)
  samp <- as.vector(sample(levels, size = nclust, replace = T))
  orig <- latam2$election
  dat <- ldply(lapply(samp, FUN = function(x){latam2[orig == x,]}), data.frame)
  
  a1<- (lm(voteshare ~ primary, dat))
  a2<- (lm(voteshare ~ primary + nparties + Decade +country, dat))
  
  return(c(coef(a1)[1:2], coef(a2)[1:4]))
}

set.seed(1234)
bs3 <- replicate(n = 2000, expr = block_bs_voteshare())
se <- apply(bs3, MARGIN = 1, FUN = sd)

sea2 <- se[1:2]
sea21 <- se[3:6]

mean(latam$voteshare)
sd(latam$voteshare)

sink("Output/Table_A2_appendix_voteshare.tex")
stargazer(vote1, vote2,
          se = list( sea2, sea21),
          dep.var.labels = "Voteshare",
          covariate.labels = c("Primary", "No. Parties", "Decade", "Constant"),
          omit = c("country+", "Constant"),
          
          omit.stat = c("f", "ser"),
          
          add.lines = list(c("Prop. of Obs. with Primaries ",rep(".23", 2)),
                           c("Voteshare Mean (SD)", rep("31.7 (16.0)", 2)),
                           c("Election FE", (c("Yes", "Yes")))),
          #  title = "Polarization - Party-Year",
          float = F)
sink()


# Figure 4 ----------------------------------------------------------------
### Equilibrium region plot

# Parameter grid
D <- seq(from = 1, to = 5, by = .01)
xm <- seq(from = -2.5, to = 2.5, by = .01)

eg <- expand.grid(D, xm) %>%
  filter(Var1 >= 2 * abs(Var2))
names(eg) <- c("D", "xm")

# One period best response functions (assuming \alpha = 12; v = 0.5, b = 2.5)

LbrRprimary <- function(Delta, xm){
  return(-0.114858 + 0.0144628 * xm + 0.053267 * Delta)
}

LbrRnoprimary <- function(Delta, xm){
  return(-0.123897 + 0.0144628 * xm + 0.053267* Delta)
}

RbrLprimary <- function(Delta, xm){
  return(-0.114858 - 0.0144628 * xm + 0.053267 * Delta)
}

RbrLnoprimary <- function(Delta, xm){
  return(-0.123897 - 0.0144628 * xm + 0.053267 * Delta)
}


op <- eg
op$LbrRprim <- mapply(LbrRprimary, Delta = op$D, xm = op$xm)
op$LbrRnoprim <- mapply(LbrRnoprimary, Delta = op$D, xm = op$xm)
op$RbrLprim <- mapply(RbrLprimary, Delta = op$D, xm = op$xm)
op$RbrLnoprim <- mapply(RbrLnoprimary, Delta = op$D, xm = op$xm)

op <- mutate(op,
             two = ifelse(LbrRprim > 0 & RbrLprim >0 , "Two Primaries", "other"),
             no = ifelse(LbrRnoprim < 0 & RbrLnoprim <0, "No Primaries", "XX"),
             R = ifelse(LbrRprim < 0 & RbrLnoprim > 0, "R Primary", "XX"),
             L = ifelse(LbrRnoprim > 0 & RbrLprim < 0, "L Primary", "XX."))

two <- filter(op, two == "Two Primaries")
two_coords <- data.frame(D = c(5, min(two$D[two$D == 2 * two$xm]), min(two$D), min(two$D[two$D == 2 * two$xm]), 5),
                         xm = c(-2.5, -1/2 * min(two$D[two$D == 2 * two$xm]), 0, 1/2 * min(two$D[two$D == 2 * two$xm]), 
                                2.5),
                         eq = "Two Primaries")

no <- filter(op, no == "No Primaries")
no_coords <- data.frame(D = c(max(no$D[no$D == 2 * no$xm]), 1, 1, max(no$D[no$D == 2 * no$xm]), max(no$D)), 
                        xm = c(-1/2 * max(no$D[no$D == 2 * no$xm]), -1/2, 1/2, 1/2 * max(no$D[no$D == 2 * no$xm]), 0),
                        eq = "No Primaries")

R <- filter(op, R == "R Primary")
R_coords <- data.frame(D = c(max(R$D[R$D == -2 * R$xm]) +.02, min(R$D[R$D == -2 * R$xm]) -.02,
                             R$D[which.max(R$xm)]),
                       xm = c(-1/2 * (max(R$D[R$D == -2 * R$xm] + .02)), -1/2 * (min(R$D[R$D == -2 * R$xm])-.02),
                              R$xm[which.max(R$xm)] +.02), 
                       eq = "R Primary")
L_coords <- data.frame(D = R_coords$D, xm = -R_coords$xm, eq = "L Primary")

coords_op <- rbind(two_coords, no_coords, R_coords, L_coords)
segs <- data.frame(x = c(two_coords$xm[2], two_coords$xm[4], no_coords$xm[1], no_coords$xm[4]),
                   xend = rep(c(two_coords$xm[3], no_coords$xm[5]), each = 2),
                   y = c(two_coords$D[2], two_coords$D[4], no_coords$D[1], no_coords$D[4]),
                   yend = rep(c(two_coords$D[3], no_coords$D[5]), each = 2),
                   eq = rep(c("Two Primaries", "No Primaries"), each = 2))


# Infinite horizon best response functions (assuming \alpha = 12; v = 0.5, b = 2.5)

R2 <- function(Delta, xm){
  return(-0.0574289 - 0.0072314 * xm + 0.0266335 * Delta)
}

L2 <- function(Delta, xm){
  return(-0.0574289 + 0.0072314 * xm + 0.0266335 * Delta)
}

R0 <- function(Delta, xm){
  return(0.124448 + 0.0181216 * xm - 0.0266335* Delta)
}

L0 <- function(Delta, xm){
  return(0.124448 - 0.0181216 * xm - 0.0266335* Delta)
}

RR <- function(Delta, xm){
  return(-0.0619485 - 0.0072314 * xm + 0.0266335 * Delta)
}

LR <- function(Delta, xm){
  return(0.113123 - 0.0181216 * xm - 0.0266335 * Delta)
}

RL <- function(Delta, xm){
  return(0.113123 + 0.0181216 * xm - 0.0266335 * Delta)
}

LL <- function(Delta, xm){
  return(-0.0619485 + 0.0072314 * xm + 0.0266335 * Delta)
}

ih <- eg
ih$L2 <- mapply(L2, Delta = ih$D, xm = ih$xm)
ih$R2 <- mapply(R2, Delta = ih$D, xm = ih$xm)
ih$L0 <- mapply(L0, Delta = ih$D, xm = ih$xm)
ih$R0 <- mapply(R0, Delta = ih$D, xm = ih$xm)
ih$LR <- mapply(LR, Delta = ih$D, xm = ih$xm)
ih$RR <- mapply(RR, Delta = ih$D, xm = ih$xm)
ih$LL <- mapply(LL, Delta = ih$D, xm = ih$xm)
ih$RL <- mapply(RL, Delta = ih$D, xm = ih$xm)

ih <- mutate(ih,
             two = 1 * (L2 > 0 & R2 > 0),
             no = 1 * (L0 > 0 & R0 > 0),
             R = 1 *(LR > 0 & RR > 0),
             L = 1 * (LL > 0 & RL > 0))

two <- filter(ih, two == 1)
two_coords <- data.frame(D = c(5, min(two$D[two$D == 2 * two$xm]), 
                               min(two$D), min(two$D[two$D == 2 * two$xm]), 5),
                         xm = c(-2.5, -1/2 * min(two$D[two$D == 2 * two$xm]), 0, 
                                1/2 * min(two$D[two$D == 2 * two$xm]), 
                                2.5),
                         eq = "Two Primaries")

no <- filter(ih, no == 1)
no_coords <- data.frame(D = c(max(no$D[no$D == 2 * no$xm]), 1, 1, max(no$D[no$D == 2 * no$xm]), max(no$D)), 
                        xm = c(-1/2 * max(no$D[no$D == 2 * no$xm]), -1/2, 1/2, 1/2 * max(no$D[no$D == 2 * no$xm]), 0),
                        eq = "No Primaries")

R <- filter(ih, R == 1)

R_coords <- data.frame(D = c(5, 
                             min(R$D[R$D == - 2* R$xm]),
                             min(R$D[R$D ==  2* R$xm]),
                             max(R$D[R$D ==  2* R$xm]),
                             5),
                       xm = c(-2.5, 
                              -1/2 * min(R$D[R$D == - 2* R$xm]),
                              1/2 * min(R$D[R$D ==  2* R$xm]),
                              1/2 * max(R$D[R$D ==  2* R$xm]),
                              max(R$xm[R$D == 5])),
                       eq = "R Primary")

L_coords <- data.frame(D = c(5, 
                             min(R$D[R$D == - 2* R$xm]),
                             min(R$D[R$D ==  2* R$xm]),
                             max(R$D[R$D ==  2* R$xm]),
                             5),
                       xm = -c(-2.5, 
                               -1/2 * min(R$D[R$D == - 2* R$xm]),
                               1/2 * min(R$D[R$D ==  2* R$xm]),
                               1/2 * max(R$D[R$D ==  2* R$xm]),
                               max(R$xm[R$D == 5])),
                       eq = "L Primary")


ih_coords <- rbind(two_coords, no_coords, R_coords, L_coords) %>% 
  mutate(mod = "Infinite Horizon")
op_coords <- mutate(coords_op, mod = "One Period")

coords <- rbind(op_coords, ih_coords)
coords$mod <- factor(coords$mod, levels = unique(coords$mod))

a <- ggplot(coords, aes(x = xm, y = D, fill = eq, col = eq)) + 
  geom_polygon(alpha = .5) + 
  theme_minimal() +
  facet_grid(.~mod) + 
  scale_fill_discrete("Equilibrium") + 
  scale_color_discrete("Equilibrium") + 
  xlab(expression(x[m])) + ylab(expression(Delta)) 

ggsave(file = "Output/Figure_4_equilibria_regions.pdf", a, width = 8, height = 4)


# Figure 5 ----------------------------------------------------------------
### Markov simulation results

mc_sim <- mutate(mc_sim, 
              twoprimaries = ifelse(twoprimaries == 2, 0, twoprimaries),
              noprimaries = ifelse(noprimaries == 2, 1, noprimaries),
              Lprimary = ifelse(Lprimary == 2, 0, Lprimary),
              Rprimary = ifelse(Rprimary == 2, 0, Rprimary),
              expprimary = twoprimaries*2+Lprimary + Rprimary,
              expprimary2 = ifelse(expprimary == 8, 0, expprimary),
              expprimary3 = ifelse(expprimary2 > 1.9, 60, 
                                   ifelse(expprimary2 == 0, 20, ntile(expprimary2, 98))))

ns<- c(1, 2, 8, 9)
brewer.pal(9, "Blues")[ns]



t2 <- filter(mc_sim, twoprimaries %in% c(1, 2) & wex >= 0 & delta == 0.05) %>%
  mutate(min_wex = wex[which.min(Delta)]) %>%
  filter(wex >= min_wex) %>%
  group_by(wex) %>%
  summarize(D = min(Delta))


segments1 <- data.frame(x = c(1),
                        y = c(2),
                        xend = c(1.3),
                        yend = c(2.6))

l2 <- lm(D ~ wex, data = t2)
segments2 <- data.frame(x = c(0.305),
                        y = coef(l2)[1]-.003 + .305 * coef(l2)[2],
                        xend =  1.242,
                        yend = coef(l2)[1] + 1.242 * coef(l2)[2]-.01)

segments3 <- data.frame(x = c(-0.013),
                        y = c(2.329),
                        xend = 1.0805,
                        yend = 2.03)

# Expected Number of Primaries
a <- filter(mc_sim, abs(wex) < (Delta/2 - .005)& wex > 0,
            delta %in% c(.05, .25, .45)) %>%
  ggplot(aes(x = wex, y = Delta, col = expprimary3)) +
  geom_point() + 
  facet_grid(~delta, labeller = label_bquote(cols = epsilon == .(delta))) +
  theme_minimal() +
  scale_colour_gradientn("Expected \nNumber \nof Primaries", 
                         colours = brewer.pal(9, "Blues"),
                         labels = c(0, .7, 1, 1.3, 2)) +
  xlab(expression(x[m]^h)) + ylab(expression(Delta)) +
  geom_segment(data = segments2, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = brewer.pal(9, "Blues")[9], lwd = 2.7)+
  geom_segment(data = segments3, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = brewer.pal(9, "Blues")[1], lwd = 2)+
  geom_segment(data = segments1, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = "gray70", lwd = 1.4)

# One primary equilibrium

c <- mutate(mc_sim, 
            oneprim = 1- noprimaries - twoprimaries) %>%
  filter(abs(wex) < Delta/2-0.01& wex > 0 & delta %in% c(.05, .25, .45)) %>%
  ggplot(aes(x = wex, y = Delta, col = oneprim)) +
  geom_point() + 
  facet_grid(.~delta, labeller = label_bquote(cols = epsilon == .(delta))) +
  theme_minimal() +
  scale_color_gradientn("Likelihood of \nOne Primary", 
                        colours = brewer.pal(9, "Blues"),
                        labels = c(0, .17, .33, .5, .67), 
                        breaks = c(0, 1/6, 2/6, 3/6, .7)) +
  xlab(expression(x[m]^h)) + ylab("") +
  geom_segment(data = segments2, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = brewer.pal(9, "Blues")[1], lwd = 2.7)+
  geom_segment(data = segments3, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = brewer.pal(9, "Blues")[1], lwd = 2)+
  geom_segment(data = segments1, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = "gray70", lwd = 1) + ylab(expression(Delta))

ga <- ggplotGrob(a)
gc <- ggplotGrob(c)
ga$widths <- gc$widths
plots <- grid.arrange(ga, gc, nrow = 2)

ggsave(filename = "Output/Figure_5_markov_chain_numerical_plots.pdf", plots, width = 8, height = 6)


# Robustness graphs: appendix ---------------------------------------------

mc_sim_appendix <- mutate(mc_sim_appendix, 
              twoprimaries = ifelse(twoprimaries == 2, 0, twoprimaries),
              noprimaries = ifelse(noprimaries == 2, 1, noprimaries),
              Lprimary = ifelse(Lprimary == 2, 0, Lprimary),
              Rprimary = ifelse(Rprimary == 2, 0, Rprimary),
              expprimary = twoprimaries*2+Lprimary + Rprimary) %>%
  mutate(expprimary2 = ifelse(expprimary == 8, 0, expprimary),
         expprimary3 = ifelse(expprimary2 > 1.9, 60, 
                             ifelse(expprimary2 == 0, 20, ntile(expprimary2, 80))))


segments1b <- data.frame(x = c(1.15),
                         y = c(2.3),
                         xend = c(1.8),
                         yend = c(3.601))

segments2b <- data.frame(x = c(0.36),
                         y = 2.765,
                         xend =  1.73,
                         yend = 3.48)

segments3b <- data.frame(x = c(-0.02),
                         y = c(2.95),
                         xend = 1.17,
                         yend = 2.35)


b <- filter(mc_sim_appendix, abs(wex) < ((Delta-0.001)/2 - .01) & Delta > 2.3 & wex > 0,
            epsilon %in% c(.05, .25, .45)) %>%
  ggplot(aes(x = wex, y = Delta, col = expprimary3)) +
  geom_point() + 
  facet_grid(~epsilon, labeller = label_bquote(cols = epsilon == .(epsilon))) +
  theme_minimal() +
  scale_colour_gradientn("Expected Number \nof Primaries", colours = brewer.pal(9, "Blues"),
                         labels = c(0, .55, 1, 1.3, 2)) +
  xlab(expression(x[m]^h)) + ylab(expression(Delta)) +
  geom_segment(data = segments2b, 
               aes(x = x, y = y, xend = xend, yend = yend), col = brewer.pal(9, "Blues")[9], lwd = 2.7)+
  geom_segment(data = segments3b, 
               aes(x = x, y = y, xend = xend, yend = yend), col = brewer.pal(9, "Blues")[1], lwd = 2)+
  geom_segment(data = segments1b, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = "gray70", lwd = 1.4)



d <- mutate(mc_sim_appendix, 
            oneprim = 1- noprimaries - twoprimaries,
            oneprim = ifelse(oneprim == 1, 0, oneprim)) %>%
  filter(abs(wex) < (Delta - 0.001)/2-0.01& wex > 0 & epsilon %in% c(.05, .25, .45) & Delta > 2.3) %>%
  ggplot(aes(x = wex, y = Delta, col = oneprim)) +
  geom_point() + 
  facet_grid(.~epsilon, labeller = label_bquote(cols = epsilon == .(epsilon))) +
  theme_minimal() +
  scale_color_gradientn("Likelihood of \nOne Primary", colours = brewer.pal(9, "Blues"),
                        labels = c(0, .17, .33, .5, .67), breaks = c(0, 1/6, 2/6, 3/6, .7)) +
  xlab(expression(x[m]^h)) + ylab("") +
  geom_segment(data = segments2b, 
               aes(x = x, y = y, xend = xend, yend = yend), col = brewer.pal(9, "Blues")[1], lwd = 2.7)+
  geom_segment(data = segments3b, 
               aes(x = x, y = y, xend = xend, yend = yend), col = brewer.pal(9, "Blues")[1], lwd = 2)+
  geom_segment(data = segments1b, 
               aes(x = x, y = y, xend = xend, yend = yend), 
               col = "gray70", lwd = 1) + ylab(expression(Delta))

gb <- ggplotGrob(b)
gd <- ggplotGrob(d)
gb$widths <- gd$widths

plots2 <- grid.arrange(gb, gd, nrow = 2)
ggsave(filename = "Output/Figure_A1_markov_chain_numerical_plots.pdf", plots2, width = 8, height = 6)

